How Many Household Formation Systems Were There in Historic Europe? A View Across 256 Regions Using Partitioning Clustering Methods

Abstract This paper reconsiders one of historical demography’s most pertinent research problems: the fiddly concept of historical household formation systems. Using a massive repository of historical census micro-data from the North Atlantic Population Project and the Mosaic project, the four markers of Hajnal’s household formation rules were operationalized for 256 regional rural populations from Catalonia in the west to central Siberia in the east, between 1700 and 1926. We then analyze these data using the Partitioning Around Medoids algorithm in order to empirically derive the “natural groups” based on the similarity and the dissimilarity of their household formation traits. Although regional differences between European household formation systems are readily identifiable, the two statistically most valid clustering solutions (k = 2; k = 4) provide a more complex picture of household formation regimes than Hajnal and his followers have been able to compile. Our finding that when regional populations cluster on similar household formation characteristics, they often come from both sides of Hajnal’s “imaginary line,” calls into question strict bipolar divisions of the continent. By and large, we show that the long-lived idea of two household formation systems in preindustrial Europe obscures considerable variability in historical family behavior, and therefore needs to be amended.


Introduction
In what has come to be seen as a totemic discourse of historical family demography, J. Hajnal (1982) compared two contrasting "household formation systems" of preindustrial societies. 1 At the core of Hajnal's proposition was a Janus-faced relationship linking age at marriage, life-cycle service, headship attainment, and household structure, which he believed accounted for most of the variation in family systems within Eurasia. When dividing the landmass into the "northwestern European simple household system" (believed to have been typical of the Scandinavian countries, the British Isles, the Low Countries, northern France, and the German-speaking lands), and the "joint household formation system" (said to have been present in the rest of the world, including eastern Europe and Asia), Hajnal stuck resolutely to vague geographic references derived from "a variety of very different conditions widely separated in time and space" (Hajnal 1982, 455). Still, he seemed to suggest that his two large-scale patterns of family systems could be conceptualized as referring to territories west and east of his famous "line" stretching from St. Petersburg to Triest (Hajnal 1965, 133).
The "two household formation systems" thesis has given rise to a voluminous body of research in which scholars have attempted to demarcate large areas of differences and similarities in demographic behaviors across Europe and beyond, and to develop various rival typologies (Smith 1981;Laslett 1983;Kertzer 1989;Duben 1990;Rettaroli 1990;Wall 1991;Smith 1993;Caftanzoglou 1994;Farago 2003;Kaser 1996;Goody 1996;Rowland 2002;Viazzo 2005;Plakans and Wetherell 2005;Dennison 2011;Head-K€ onig and Pozsgai 2012;Szołtysek 2008Szołtysek , 2015also Mitterauer 1999;Ruggles 2009). By the end of the 1990s, it was increasingly suspected that the empirical regularity of Hajnal's markers is subject to many gradations, and that the geographic distribution of these markers may not allow for clear-cut divisions of the continent (e.g., Wall 1991;Kaser 1996;Rowland 2002; also Dennison and Ogilvie 2014;Micheli 2018). The sharpest criticism of Hajnal's demarcation of European household formation systems came, however, from the Italian and the Iberian peninsula, where an unexpected degree of regional variation in the relationship between household structure, marriage, and the institution of service has been reported, indicating that Hajnal's thesis is applicable neither in Italy, nor in Spain or Portugal (Da Molin 1990;Barbagli 1991;Viazzo 2005; also Kertzer and Hogan 1991;Manfredini 2016;Kertzer and Brettell 1987;Rowland 2002, esp. 66-69).
More than 20 years after these seeds of doubts were first planted, the most crucial implications of this critique have yet to be tackled, and the allimportant classification questions remain unresolved. Although scholars have been increasingly suspicious of Hajnal's two household formation systems theory as being too coarse-grained, most of the empirical counter-evidence remains regionally-bounded. This begs a crucial question whether the claims made about Italy or Portugal may have wider European relevance, or they can be dismissed as a mere peculiarity of the Mediterranean. However, until very recently, a rigorous examination of Hajnal's thesis across broader European terrain was impossible because it would require harmonized large-scale data infrastructures and methodologies that did not yet exist, or that were not known to family historians. While some local and small-scale regional studies attempted to test the markers of Hajnal's typology beyond Italy, most of these studies proved ineffective at improving the clarity of the cartography of European household formation patterns (e.g., Duben 1990;Caftanzoglou 1994;Smith 1993;Schlumbohm 1996). For many areas of Europe, hardly any data were available, which made systematic comparisons of household formation patterns and exploration of their spatial contours difficult to undertake (e.g., Szołtysek 2012). Hajnal treated his "household formation rules" as a complex system of intertwined elements that should be analyzed conjointly, but gaps in the historical evidence have prompted scholars to study one of the "Hajnalian" markers while neglecting the others based on their own predilections or data availability (e.g., Da Molin 1990;Barbagli 1991;Caftanzoglou 1994; see also Dennison and Ogilvie 2014;cf. Smith 1993).
In this paper, we set out to overcome these difficulties by reconsidering Hajnal's original proposition using unprecedented historical data infrastructures from both sides of his "imaginary line" and formalized tools of data mining and pattern recognition. Our primary goal is exploratory, and involves addressing the following questions: Do European historical populations form any natural groupings based on how similar or dissimilar they are with regards to Hajnal's household formation markers of marriage, service, headship, and household structure patterns; and, if so, how many such groupings can plausibly be identified? 2 Is Hajnal's dual typology of household formation rules valid for historic Europe, or are there other taxonomies in which the tradeoff between information entropy and precision is better? Does the spatial distribution of these groupings support the claim that the European continent can be divided into two distinct household formation regimes along the east-west axis? Finally, are markers of household formation systems invariably associated with each other in different areas, as Hajnal's hypothesis seems to imply?
Chasing answers to these questions is very timely as opportunities now exist to rid the classification of historical household formation patterns of its traditionally speculative nature. The recent information revolution in historical population studies (Ruggles et al. 2011;Ruggles 2012;Szołtysek and Gruber 2016) has made historical demographic data more broadly available. Rigorous data harmonization schemes allow historians to effectively measure all Hajnalian household formation markers across multiple locations, and over time. Finally, the last two decades have witnessed enormous developments in data mining methodologies, whereby techniques known as unsupervised machine learning or cluster analysis (e.g., Everitt, Landau, and Leese 2001;Hastie, Tibshirani, and Friedman 2009) allow scholars to derive optimal groupings in multidimensional data with high efficiency but little cost.
We expand on the existing literature in four major ways. First, we combine data from the North Atlantic Population Project and the Mosaic project to create an unprecedented family history database covering 256 rural populations from Catalonia to Siberia and from Tromsø to Albania. Second, by explicitly acknowledging Hajnal's formulation of household formation rules as a set of four interrelated axioms, we identify Hajnal's critical markers and operationalize them by means of harmonized measures, which are then analyzed conjointly across hundreds of Eurasian regions. Third, instead of relying on the ad hoc deductive typologies that have dominated previous studies, we employ formal methods of automatic (unsupervised) pattern recognition (Partitioning Around Medoids clustering algorithm; henceforth PAM) to discern the optimal natural groupings among our populations based on the similarities and the dissimilarities in their household formation markers. Finally, we push the frontiers of previous research by exploring the relationships between the markers of Hajnal's household formation systems among the obtained partitions.
Our results show that although regional differences between European household formation systems are readily identifiable, the two most statistically valid groupings (two-and four-cluster solutions) provide a more complex picture of household formation regimes than Hajnal or other researchers that followed him have been able to assemble. Our final cluster-based taxonomy suggests that there are four distinct patterns of household formation, two of which are impossible to classify within the Hajnalian framework. Hajnal's model is also called into question by the evidence we found that when regional populations cluster on similar household formation characteristics, they often come from both sides of Hajnal's "imaginary line." Last but not least, we show that various "scalar types" of the associations between the four elements of household formation systems are possible, and that the systemic relationship between marriage age, lifecycle service, headship attainment, and household structure can vary between clusters.
The remainder of this paper is divided into four main parts. In section 2, we present our data and methodology. We also provide a rationale for our chosen clustering strategy, and an assessment of the optimal clustering solutions that are applied to our data. In section 3, we summarize our results by presenting the statistical and the spatial distribution of the clusters of household formation markers we have identified, and describing their characteristics. In section 4, we test the sensitivity of our clustering results to different data perturbations. In closing, we discuss the possible implications and limitations of our study, and suggest some research agendas for the future.

Data
We use historical census micro-data from the harmonized and geo-referenced databases of the North Atlantic Population Project (NAPP) and the Mosaic project, which have already been discussed in detail on several occasions before (see Szołtysek and Gruber 2016;Szołtysek et al. 2017;Szołtysek and Poniat 2018; see also Online Appendix). 3 A novelty of the current version of the dataset is that it includes historical census micro-data from Siberia that cover a large share of the indigenous peoples of Russia's circumpolar north (Anderson 2011), thus making a first ever foray into northwest Asia (see Figure 1A,B). Since Hajnal related his theory of household formation systems to rural populations (1982,(450)(451), a rural dataset from the overall NAPP and Mosaic collections has been derived. 4 In total, we use in this study Mosaic samples containing 708,564 individuals living in 126,349 family households in 106 geographic areas. The NAPP samples include data for 150 additional regions from five national censuses that cover more than 11 million individuals in nearly 2.5 million domestic groups. 5 The combined 256 regional populations span the period from 1700 to 1926/27. Of these historic population samples, 22.7 per cent predate 1800, 17.6 per cent are from between 1800 and 1850, while 59.8 per cent are from after the mid-19th century.
Our approach is situated at the meso (regional) level of comparative analysis (see Hajnal 1965Hajnal , 1982Laslett 1983;Wall 1995;Manfredini 2016). 6 Accordingly, all of our focal variables are defined as regional rates (population means or age-specific rates). Furthermore, the regional data are treated as pooled time cross-sections based on the tacit assumption that the family behaviors they pertain to represent "deep" cultural layers that move slowly over time (e.g., Reher 1998;also Sch€ urer et al. 2018). While this approach has been promoted by the doyens of comparative family history, including Hajnal (see more in the Discussion and Conclusion section), it is further justified in our case because the overwhelming majority of our regional populations (240 out of 256) represent pre-transitional demographic regimes (see Szołtysek et al. 2019). Hajnal (1982) suggested that there was a functional relationship linking marriage, premarital life-cycle service, and entry into headship; and argued that one or more of these powerful variables led to the emergence of simple household systems in north-western Europe, and to joint (complex/grand) family systems in other areas. Thus, he seemed to be arguing that the European household formation systems included both the patterns themselves, as well as the modes of behavior that created those patterns. Accordingly, our operationalisation of his household formation markers includes a set of four variables, all of which we computed from the harmonized historical census micro-data used in this paper: female age at marriage, the incidence of female premarital service, the relationship between marriage and headship attainment, and the share of nuclear households (also Smith 1993, 330;Smith 1981, 600). 7 Female age at marriage: The female age at marriage was calculated using Hajnal's formula for SMAM, which relies on the age-specific distribution of people by marital status between ages 15 and 54.

Operationalization of Hajnal's Household Formation Markers
Female service: the incidence of pre-marital life-cycle service is captured with a proportion of unmarried female servants among unmarried women in the age group that is determined by the value of the female SMAM in each of populations under study. 8 For regions where female SMAM is equal to 22.5 years or lower (47 regions), the proportion of servants in the 15-19 age category is considered when computing this variable. For the populations in which the age at first marriage ranged from 22.5 to 27.5 (157 regions), service is measured in the 20-24 age group. In all of the remaining cases (female SMAM greater than 27.5-52 regions), the proportion of servants was derived for the population between ages 25 and 29. This variable provides a better approximation of the phenomenon Hajnal was referring to than more traditional measures, such as the percentage of servants in the population.
Proportion of simple families: This is the regional share of simple (nuclear) households, as derived by applying the Hammel-Laslett classification (Hammel and Laslett 1974, 92) to our dataset, in which nuclear households are defined as consisting of one conjugal family unit only. While we acknowledge that this measure may not be sufficiently meaningful for more sensitive analyses of household systems (Berkner 1972;Ruggles 2012), we believe it serves the present purpose because it distinguishes between simple and multiplefamily households, and thus adheres closely to Hajnal's original formulation.
Marriage-headship nexus: Hajnal proposed to measure the relationship between marriage and headship attainment by comparing the age-specific curves of Spatial data distribution by type of data collection and major institutional and socioeconomic groupings of respective populations. Notes: each point on the map represents one Mosaic/NAPP regional population as defined in the text. Seven bigger territorial groupings on the right-side panel of the figure followed major institutional and socioeconomic distinction across historic Europe. "Great Britain": England, Wales, and Scotland; "Scandinavia": Danish, Swedish, and Norwegian data, as well as Iceland; "Germany": German dominated areas other than the Habsburg territories; "West": areas west and south-west of Germany; "Habsburg": Austrian, Hungarian, Croatian, as well as Slovakian data; "East": east-central and eastern Europe, including the former Polish-Lithuanian Commonwealth and Russia (including Siberian territories geographically in Asia); "Balkans": areas south and/or east of Croatia and Hungary. Source: Mosaic/NAPP data. See Online Appendix for the list of all Mosaic/NAPP datasets. entry into marriage and into headship (Hajnal 1982, 463). Because attempting such a comparison for all 256 regions would be unproductive and burdensome, we derive the synthetic measure called the Cumulative Marriage Headship Difference (CMHD). Given that the marital age among the men in our data generally ranges from the twenties to the mid-thirtiesand that for the majority of our populations, male headship rates converged toward a plateau among men in their fortieswe decided to consider all age groups between ages 20 and 40 for the computation of this measure. Accordingly, the CMHD is an outcome of the subtraction of the areas under the two curves, with one reflecting the age-specific proportion of ever-married individuals, and the other reflecting the respective proportion of ever-married heads. The area under the proportion of ever-married curve can be understood as the mean number of years lived by the member of a synthetic cohort as ever-married, and the area under the proportion of ever-married heads curve as the mean number of years lived as the ever-married household head. Accordingly, the difference between those two values is the average time between getting married and attaining the headship in the 20-40 age classes. The areas under the curves (definite integrals) are approximated using the trapezoidal rule. 9 Raw Data Spatial Distribution Figure 2 presents a descriptive mapping of the four variables selected for clustering. When looking at this figure, we can clearly see that the bipolar division of the continent cannot be generalized. The figure shows, for example, that the female SMAM remains highest in the broadly defined north-western and western parts of Europe. But the true axis with the highest female ages (27-31 years) runs somewhat vertically through Norway and southern Sweden, and then through much of the German-speaking lands, down to Switzerland and Austria. While female marital ages are lower on both sides of this apparently broad zone, they are much lower only to the east and southeast of it. However, the data for the eastern zone are far from uniform. Whereas the lowest-low women's ages run crosswise from southern Belarus, central Ukraine, and the eastern Balkans toward Albania; the values in the remaining part of the eastern zone are much less extreme. Of these areas, the medium low pattern of East-Central Europe and much of Russia (20-24) seem to be most generalized. But even within the Russian territory, there are regions where the marriage patterns are much more similar to those observed in the west (these populations come from both ends of our chronological span).
When we look at the female service variable, we can detect three general patterns. First, a very high service incidence area can be observed through parts of Scandinavia (especially in Norway and Denmark, though to a much lesser extent in Sweden), with some additional hot spots dispersed over Latvia (Kurland), central Poland, north-western Germany and the Netherlands. The rest of our data are split between a large area with more moderate values that covers much of England, Wales, and Scotland, parts of France and the Low Countries, and certain parts of East-Central Europe; and a long but discontinuous area with very low values at the southern and eastern edges of the map, particularly in the Balkans and to the east of present-day Poland. An exception to this general pattern is the coincidence of comparably low service areas in Switzerland, parts of France, and Catalonia.
While the distribution of the proportion of nuclear families follows the distribution of the female SMAM to some extent, the overall pattern is spottier. Again, we can see a vertical belt of a high incidence of household simplicity running across the data, albeit with more pronounced intrusions into central Europe (especially Poland) and even into the Balkans (Romania), where it appears that the populations had the same high propensity for household simplification as the populations in many Swedish, Danish, Belgian, and German regions, as well as in Iceland. Again, to the east and the south of this area, we can see a substantial decline in the share of nuclear households, particularly in the Baltic area, Belarus, Albania, and Serbia. While the pattern in Russia generally converges with this trend, it is still quite diverse, as Russia's Siberian populations seem to have had more streamlined household structures. We also see that while the populations of England and Wales were not at the forefront of structural simplicity, they still tended to favor household nuclearity, and to a greater extent than the Scottish populations in the north. Finally, we note that populations in which the share of nuclear households was low medium (0.47-0.64) are observed across the various macro-regions of Europe, including in France, Scotland, German-speaking areas, Ukraine, and even Siberia.
When we consider the CMHD, we see that in nearly all of the north-western and western populations, the average time between getting married and attaining headship was very short. Nearly everywhere else, the relationship between marriage and headship was much weaker. While we can see that the CMHD patterns in Poland, Hungary, and even part of Ukraine are similar to those in France, Catalonia, and north-western Germany, we note that the cumulative gap between marriage and headship increases substantially to the east of the Bug river and in the Balkans. However, these patterns still diverge from the Hajnalian divide, as populations with marriage-headship intensities similar to those in north-western and western Europe are shown to be present in Romania/ Wallachia and eastern Siberia.

Clustering Algorithm
Taken together, the observations presented above suggest that the distribution of the focal variables intermingled to a large degree in terms of both their quantities and their spatial configurations. A major conceptual challenge that arises in this context is that various combinations could be imagined and developed along these four dimensions to describe the historical household formation systems across our 256 populations. Cluster analysis, which allows for the inclusion of multiple variables as a source of configuration definition, seems particularly useful in this situation. Through cluster analysis, we can group our populations into classes in which the degree of similarity is high between members of the same class, and is low between different clusters (Everitt, Landau, and Leese 2001)without any a priori expectation about the data's underlying structure (Han, Kamber, and Pei 2012, 444;Ahlquist and Breunig 2012). 10 Accordingly, the 256 populations were subjected to the robust Partitioning Around Medoids (PAM) clustering algorithm (Kaufman and Rousseeuw 1990;Kassambara 2017, 48-56;also Jin and Han 2017). 11 There are three main reasons why we chose to use PAM. First, the nonhierarchical PAM method was appropriate because in light of our domain knowledge, there was no reason why a hierarchy of clusters should be imposed on our data. PAM assigns objects to mutually exclusive clusters based on certain predefined criteria after the number of clusters to be formed is specified prior run. The end result is a number of disjointed, nonoverlapping groupings that together make up the full dataset. Second, the mechanics of a partitional clustering algorithm such as PAM allowed us to minimize the degree of subjective judgment that is often implicit in the interpretation of the hierarchy of nested clusters generated by hierarchical methods (Zambelli 2016). Third, PAM outperforms other unsupervised partitional clustering algorithms, such as k-means (Jain 2010), that are highly sensitive to outliers (Han, Kamber, and Pei 2012, 454;Kassambara 2017, 48), by choosing medoids rather than means as the cluster centroids; where a medoid is the data object of the cluster that is most centrally located (Ayramo and Karkkainen 2006). 12 Furthermore, partitional clustering algorithms like PAM have advantages in applications involving datasets for which the construction of a dendrogram may be computationally and visually prohibitive, such as in hierarchical clustering (Jain, Murty, and Flynn 1999, 278). 13 Typically, PAM (like k-means) is run independently for different values of k (the user-specified number of clusters; 10 in our case), and the partition that appears to be the most robust based on a range of complimentary statistical test for the optimal clustering solutions is chosen. Our own determination of the optimal number of clusters was based on three widely used methods: average silhouette widths, total withincluster sums of squares (so-called Elbow plot), and the Gap statistics (Kaufman and Rousseeuw 1990;Tibshirani, Walther, and Hastie 2001;Halkidi, Vazirgiannis, and Hennig 2015;Kassambara 2017, 128-137). 14 These validation measures converged toward a limited number of preferable clustering solutions. Both the Elbow and the silhouette methods indicate that using two clusters is the optimal choice for our data, but that using k ¼ 4 is the second-best choice. 15 However, the Gap statistic, which is considered to be the most robust and objective method (Tibshirani, Walther, and Hastie 2001;Hastie, Tibshirani, and Friedman 2009, 519-520;Kassambara 2017, 135), strongly indicates that using a four-cluster solution is the optimal choice for our data.
Independent of the mathematical formalism argument, the suggestion that either the two-cluster or the four-cluster structure is necessary and sufficient to explain the observed patterns of household formation variation in our data seems intuitively appealing (see (Hennig 2015). These structures enable us to compare two seemingly equivalent categorization approaches: the approach that is most similar to Hajnal's original bipolar partition of the continent; and a more nuanced approach that goes beyond this division, and thus allows us to amend Hajnal's original proposition with revisionist stances that emphasize the regional peculiarities of family organization (see also the Discussion and Conclusion section).

Clustering Results
In this section, we proceed with visualizing the clustering results and mapping them onto the cartograms. Their derivative groupings are further examined with respect to their internal characteristics and cohesion. It should be noted that the maps shown in this paper represent only one way that the outcomes of statistical clustering could be presented. The k-medoids algorithm does not take into account the spatial proximity or contiguity. While we are aware of formal spatial statistics that would allow us to identify nonrandom spatial clusters (Anselin 1988), given the aim of our study, we have chosen to use procedures that ensure that our results can be compared with the traditional, spatiallyinsensitive classifications of earlier scholarship.

Partition of the Hajnalian Variables on Two Cluster Solutions
The partitioning around two medoids (k ¼ 2) results in two clusters of substantially different sizes. Cluster 1 (tentatively called "eastern") contains 45 regions, while Cluster 2 (tentatively called "western") includes 211 regions. In order to visualize this cluster partition statistically, dimensionality reduction through the Principal Component Analysis (PCA; Everitt, Landau, and Leese 2001, 29-32) was applied. The two new variables derived were then used as the basis for the plot according to the first two principal components that explain the majority of the variance that is expressed by all the four original input variables (Figure 3). 16 Figure 3 shows that the obtained partition is pretty cleanly separated. There is no appreciable overlap between the clusters, and the clusters are roughly similar in shape. This geometry is reflected in the summary statistical characteristics of the two groupings (Table 1; Figure 4), as the median values of the four Hajnalian markers diverge sharply between the two partitions. A possible interpretation of the "western" cluster is that it exemplifies the Hajnalian "north-western pattern," which is characterized by the average female age at marriage above 26 years; a remarkable incidence of     female life-cycle service (somewhat 13 times higher than it is in the "eastern" cluster); the very compressed time between entry into marriage and entry into headship, and the domination of nuclear households. Accordingly, the medoid of the "eastern" cluster is in many respects the mirror image of the previous cluster, and is very much in line with the "joint" system originally envisaged by Hajnal. In this cluster, female marriage occurs much earlier, the incidence of premarital service among women is very low, marriage-headship relationship is loose, and an average proportion of nuclear families is below 50 per cent.
The parallel coordinate plot, in which the difference from the mean for each cluster medoid is provided for all of the constitutive variables in the partitioning ( Figure 5), further supports this observation. It shows that for all four dimensions, the two clustered populations differ considerably around the key four characteristics, albeit most strongly with regard to the marriage-headship nexus. Taken together, these findings suggest that the partitioning on two solutions provides a powerful distinction between different demographic regimes, and hence offers a high degree of consistency with the Hajnalian hypothesis of two household formation systems.
However, a clustering solution that seems convincing on purely mathematical grounds may not look as attractive when it is mapped onto the geographic coordinates of the member populations ( Figure 6). At first glance, it appears that the bipolar division of the continent is a very accurate reflection of reality, as most of the observations from Cluster 2 are located in the west and the north-west, as well as on the Scandinavian rim of our data; and that the majority of the data points from Cluster 1 are located mainly in the Balkans and to the east of the contemporary border between Poland and former Soviet republics. However, there are three major caveats to this apparently dichotomous rendering of the geography. First, Cluster 2 extends well beyond the north-western "core" areas to which it is supposedly confined, spanning from Iceland to the Romanian (Wallachian) shores of the Black Sea, and from Catalonia to central Poland. This pattern indicates that in some parts of central and east-central Europe (e.g., in Poland, Romania, Slovakia, Hungary, and parts of western Ukraine) where the "eastern" household formation pattern was assumed to dominate in Hajnal's classification and in the literature it has inspired, the household formation markers were actually much more similar to those in England, Germany, and Scandinavia than to those in the Balkans and Russia.
Other important caveats are that the "eastern" cluster is split between two largely disjointed areas in south-eastern Europe (without Romania) and east of Poland; and that some of the easternmost populations in our datarepresented by the Siberian indigenous tribes of Tungus and Yakuts on the cartogramare also members of the dominant "western" cluster. Therefore, while a general east-west contrast remains, it is obvious that an absolute east-west division among the members of the two clusters is difficult to impose. Breaking up the macro-regional groupings within our data (see Figure 1B for definitions) by their cluster membership further supports this observation, as it shows that one-quarter of all data points from the "east" and more than one-third of all data points from the "Balkans" are allocated to the dominant "western" cluster. It should also be noted that the above observations are generally robust to different data perturbations (see section Sensitivity Tests). Thus, when measured against Hajnal's thesis, these results suggest that even if household formation dualism seems to be reified in statistical terms across our data, it is much harder to demonstrate in strictly geographical terms. Therefore, Hajnal's own prediction that north-western European household formation systems are likely to be found among populations outside of north-western Europe (see Hajnal 1982, 450) appears to be confirmed.

Partition of the Hajnalian Variables on the Four-Cluster Solution
As might be expected, the four-cluster solution offers a much more nuanced view of the distribution of household formation markers, suggesting that Hajnal's dichotomous view masks important variations on both sides of his "imaginary line." Introducing additional clusters 17 results in a stronger differentiation of a formerly monolithic two-tiered pattern, even though the alterations introduced by the new partitioning approach have different effects in different portions of our data (see Figure 7).
The "east" cluster from the two-class partition remains partly intact. Its prime areas of concentration continue to be in the Balkans, the territories of present-day Ukraine, Belarus, the Baltic countries, and the European part of Russia; and with no parallels whatsoever in the areas to the west of the Bug river and to the north of Belgrade. However, unlike in the two-cluster structure, its coverage of the easternmost part of our data becomes much more idiosyncratic, as populations to the east of the Urals break sharply from the patterns that dominate the European part of Russia. This result suggests that while the "joint family"-like household formation pattern dominates in the east, it does not apply universally across eastern Europe and northwest Asia, and that its geographical reach is delineated by a double "fault line": one that demarcates the European east from the west, and another that demarcates the western areas of this eastern zone from the Siberian areas to the east. This observation is further supported by the finding that the "east" region was the only one of the seven major institutional and socioeconomic macro-regions in which the populations are represented across all four clusters, as we show here.
Even more pronounced changes occurred within the earlier "western" cluster. Figure 7 shows that whereas there was previously a largely uniform bloc stretching from Iceland to Romania, and from Catalonia to Tromsø and Murmansk, this bloc has now been split into three separate groups. England, Wales, and Scotland stand alone as the most uniform spatial conglomerate, but the pattern encapsulated by the populations of these countries is by no means unique. It extends beyond the British Isles to cover much of France; the Low Countries; the northern, western, and southern Mosaic areas of Germany; Switzerland and Austria; andafter passing through south-western Polandalmost all of Sweden in 1880 (now Cluster 3). Two broad regions in central Siberia with indigenous peoples are the "strange bedfellows" in this otherwise westernoriented grouping. In fact, Cluster 3 is now the only cluster in which the members are spread over six out of the seven major institutional and socioeconomic macro-regions that encompass our data collection. 18 Within the northern-central rim that was previously occupied by the "western" cluster, a new partition -Cluster 4has been identified. This cluster appears to be the most cohesive in spatial terms (except for Iceland). It divides Scandinavia by covering Norway and Denmark, but not Sweden (except for the south-eastern coast of the country), and then crosses the Baltic Sea to make forays into a number of areas in northern and central Germany, and even further south-east, toward northern and central Poland. Still, in our dataset, the spatial range of this cluster never seems to reach beyond the North European Plain in the south, or beyond the Vistula river in the east.
A final deviation from the two-cluster solution involves the delineation of the new Cluster 2. Strong regional boundedness could have been ascribed to this grouping, which primarily covers southern Poland, Slovakia, Hungary, Romania, and western Ukraine. However, the cluster also covers populations located quite far away from its centroid (located in Podolia in the Ukraine), including Catalonians and the majority of the Siberian populations in the east. Our finding that there are clear affinities in household formation rules among societies so diagonally dispersed represents another obvious challenge to the rigid spatial grouping formulated by Hajnal.
Having discussed the spatial contours of the fourcluster partition, we display in Figure 8 the size and the shape of the clusters, as well as their relative position (as before, a dimensionality reduction algorithm via PCA was used). 19 Clusters 2-4 are quite tight, but Cluster 1 is somewhat more diffused. An inspection of the cluster plot generally reports a satisfactory conclusion, as a majority of the clusters are non-overlapping and do not have joint intersections. However, Cluster 3, which only minimally overlaps with Cluster 2, overlaps more with Cluster 4, what could indicate a potentially ambiguous relationship. Viewed together, Figures 7 and 8 suggest that despite their formal distinctiveness, the populations in Cluster 1 and Cluster 2, as well as the populations in Cluster 2 and Cluster 4, can be spatially adjacent (like in central Poland) or even contiguous (e.g., in central Ukraine).
The patterns that underlie this geometry can be further explained by looking at which variables (and in which direction) are driving the four-cluster structure. Compared to the dual partition, Figure 9 shows more complex patterns. As expected, the attribute values are the most divergent for clusters 1 and 4, confirming their strongly antipodal character. Clusters 3 and 4 tend to diverge most on the importance of female service, but they otherwise converge on female marriage andalbeit to a lesser degreeon marriage-headship nexus, which explains the partial overlap revealed in Figure 8 above. Another interesting pattern that resurfaces is that while the values for Cluster 2 are strongly diametric to the values for Cluster 4 on female premarital-service and marriage, they are more aligned on the incidence of nuclear households. However, the relative closeness of clusters 1 and 2 observed for the marriage and service variables largely disappears when the headship and  household patterns of the two partitions are compared (also Figure 10 and Table 2).

The Amended Taxonomy of Household Formation Systems
Taken together, these characteristics combine to produce the amended taxonomy of household formation systems derived from the four-cluster structure: CLUSTER 4: This cluster represents an extreme example of the Hajnalian north-western European simple households system. It is characterized by very late female marriage (a median age above 27); a very high incidence of premarital service among unmarried women (close to 60 per cent); strong neolocality, as encapsulated by an almost absolute marriage-headship nexus, and the highest average proportion of nuclear family households across all partitions (median of 76 per cent). All of the Danish populations from the 1787 census (except one) and all but one of the Norwegian populations from the 1801 census are assigned to this cluster, together with the populations of 18th-century Iceland and three Polish regions, but only two English and two Swedish regional populations. It therefore seems ironic that the cluster does not include the countries on which the very idea of a "unique" north-western European family pattern was habitually modeled; i.e., England and the Netherlands (see Laslett 1977;Hendrickx 2005; and, for a similar perspective on this point, Dennison and Ogilvie 2014). Cluster 4 consists primarily of populations surveyed before the mid-19th century, but is otherwise not particularly coherent chronologically. Most of its members were only minimally exposed to the potential effects of larger industrialization processes.
CLUSTER 1: Its characteristics are very much akin to the Hajnalian joint (complex) household formation category. This cluster stands out as having the lowest age at marriage for women, with a median SMAM value of 20.8 years; the lowest median proportion of  premarital service (albeit with some prominent outliers); a very loose, if any, relationship between entry into marriage and entry into headship; a very low incidence of nuclear households (around 40 per cent on average, but with more than 50 per cent of observations well below that threshold, and never exceeding 56 per cent). Despite its stark statistical distinctiveness, Cluster 1 actually has weak geographical (and chronological) coherence, as its two branchesi.e., the Balkan branch and the eastern European branchare divided spatially by the presence of Cluster 2, which wedges in between them. Yet despite this incoherence, all of the member populations of this cluster had not yet been affected by the start of the industrial revolution.
Statistically speaking, these two clusters encapsulate the opposing features that Hajnal seemed to have in mind when contrasting household formation systems. Indeed, they include some of the populations Hajnal referred to (1982,456,467) when envisaging this dichotomy, such as Denmark and Russia. However, in addition to these structural-demographic antipodes, two other distinct patterns can be identified.
CLUSTER 3: Its core areas spread over Great Britain, France, the Netherlands, Belgium, much of Sweden, some of the German-speaking areas (including Austria), Switzerland, and south-western Poland. Compared to Cluster 4, this cluster has a similar pattern of late female marriage (median age of 26.5), and an only slightly less pronounced tendency toward neolocality. It might have been seen as forming a bigger "family" together with Cluster 4 (see above) had it not departed from the latter in another two dimensions. While Cluster 3 still has a much higher rate of service than Clusters 1 and 2, the share of servants among unmarried women in Cluster 3 is nearly half that of Cluster 4. The share of nuclear households is also much lower in Cluster 3 than in Cluster 4, with the lower 50 per cent of the former's members falling well below the lowest quartile of Cluster 4. This difference can be explained in part by the observation that Cluster 3 comprises a number of populations that are known to have leaned toward stem family organization, and were thus treated rather ambiguously in Hajnal's hypothesis (see Engelen and Wolf 2005, 16-18), such as populations in north-western Germany (Westphalia) and southern France. 20 In conclusion, we are inclined to see Cluster 3 as representing a mixture of populations, including populations that deviated from some of the markers of Cluster 4, but that otherwise had the same essential features; 21 as well as populations that had qualitatively different modes of household formation, such as those based on stem family principles. CLUSTER 2: This cluster represents a vivid manifestation of how the apparently dichotomous household formation markers suggested by Hajnal could combine and display enough statistical coherence to be recognized as a "natural grouping" by the clustering algorithm. Two features of this group of populations appear to be in line with Hajnal's understanding of the characteristics of joint families and of eastern household formation rules. One of these characteristics is an early female age at marriage (median of 21.4 years). While this may suggest the populations of Cluster 2 closely resemble the populations of Cluster 1, this similarity is far from absolute, given that the lowest quartile of the marriage age distribution of the former falls well below the overall spread of the data from Cluster 2, as shown in Figure 10 above. The low proportion of servants among unmarried women is another "joint family"-like trait that is also shared with the populations of Cluster 1. Meanwhile, this grouping differs from the other clusters in that it combines those characteristics with features that either fit the Hajnalian "north-western European" standard, or even encapsulate the standard's extremities. The CMHD pattern of Cluster 2 exemplifies the former tendency. Whereas its median value is roughly halfway between that of clusters 1 and 4, the lower 50 per cent of this attribute's distribution overlaps with the absolute range of Cluster 3, while the lowest quartile overlaps with the spread within Cluster 4. Given that the CMHD patterns we have described imply that neolocal household formation was popular, if not prevalent, among these populations, it is not surprising that as in groupings 3 and 4, the share of nuclear households is also very high in this cluster. It should, however, be noted that the shares of nuclear households in the upper quartile fractions of the Cluster 2 generally override the highest values ever encountered in these other clusters (values above 82 per cent), which suggests an extreme simplification of the household structure.

Sensitivity Tests: Cluster Validation
To test the robustness of our taxonomy, the above partitions of our dataset were evaluated using selected cluster validity methods (Theodoridis and Koutroubas 1999;Halkidi, Vazirgiannis, and Hennig 2015). Because there is in our case no "true" classification of household formation systems that is assumed to exist, and against which our results might be validated, the fit between the structures imposed by the clustering algorithm and the data is tested using the data alone (Jain 2010, 657).
One way to assess whether a cluster represents true structures is to see how sensitive it is to perturbations in the input data (Hennig 2007). In the first step, we tested the sensitivity of the cluster structure (both k ¼ 2 and k ¼ 4) to the elimination of variables (attributes) using four different stability measures, such as the average proportion of non-overlap (APN), the average distance (AD), the average distance between means (ADM), and the figure of merit (FOM) (Kassambara 2017, 152). 22 All of these methods are based on the cross-classification table of the original clustering, with the clustering based on the removal of one column (attribute). While the k ¼ 2 approach scores better on APN and ADM, k ¼ 4 outperforms the former approach on AD and FOM. In all cases, the differences between the two-cluster and the four-cluster solutions are small, which suggests that the two partitions are similarly stable.
In the next step, we tested the sensitivity of the final cluster structure when subsets of observations (rows), rather than attributes, are changed, using the Jaccard coefficient, which distinguishes between stable and spurious clusters (R library fpc, Hennig 2018). 23 The obtained results suggest that the cluster stability of the four-structure partition is high. Even the worst-performing Cluster 2 always yields a mean Jaccard coefficient above 0.76, where 0.75 is considered a benchmark for a valid, stable cluster (Hennig 2007). Other clusters have much higher levels of similarity (generally, well above 0.85, and most often above 0.90).
In addition, we tested the robustness of our partitions by removing Great Britain from our sample; and, second, by replacing the female variables for pre-marital service and marriage with their male equivalents. Given that the English (1851) and the Scottish (1881) rural data from our sample represent populations from territories where industrialization was already well advanced, it is possible to question whether these populations truly reflect the "traditional" social structures that Hajnal had in mind (e.g., Gritt 2000;Sch€ urer et al. 2018). Accordingly, re-running the clustering algorithm with the malecentred variables allows checking how sensitive our results might be to different developmental trajectories of male and female service and marriage (e.g., Hinde 1985).
However, both tests generally produced results very similar to those of the base models. The partitioning without England, Wales, and Scotland yielded almost identical patterns, especially for the four-cluster partition (cartograms available on request). Re-running the clustering algorithm with male variables also revealed a high degree of overlap with the base findings. The Adjusted Mutual Information index (AMI; Vinh, Epps, and Bailey 2010) is 0.80 for k ¼ 2. For k ¼ 4 the AMI is lower, but still reasonable (0.48), mainly due to the fact that in the new partitioning most of England (unlike Scotland and Wales) has been allocated to Cluster 2 from the base findings (cartograms available on request). Yet in all other respects, the basic spatial contours of this perturbed partition were very similar to those of the base model.
Altogether, the results from these exercises offer a second level of confirmation to the basic inferences drawn from the cluster analysis.

Discussion and Conclusion
Carrying on from where Mediterranean scholars stopped some 25 years ago, this paper has reconsidered the concept of two household formation systems formulated by J. Hajnal using a hitherto unprecedented collection of historical census microdata and data mining techniques previously not applied by family historians. The empirical findings of this article can be summarized as follows. Readily identifiable regional differences (statistical clusters) between household formation systems have indeed existed across Europe, but widening the observation of Europe's household formation patterns to 256 regional populations revealed considerably more variation than Hajnal and his followers had predicted. Although Hajnal's dual partition concept has received considerable statistical support, we showed that the geographic framework of this division was not robust enough to absorb substantial counterevidence without changing its main thrust; and, above all, that it raised the comparative discussion of European household formation systems to such a high level of aggregation that on-the-ground diversity was never properly acknowledged. Most crucially, based on careful optimization criteria, we found evidence that the fourcluster division of household formation regimes offers a far more reasonable way of capturing the variation across our data than the dual partition model, providing the best tradeoff between taxonomic distinctions and the acknowledgement of the local peculiarities of family organization. Mathematical formalism aside, the four-group clustering also has an advantage of reconciling the dissonant and dispersed perspectives of a number of scholars who previously questioned the undifferentiated nature of Hajnal's model (e.g., Kertzer 1991a;Wall 1991;Farago 2003;Szołtysek 2008Szołtysek , 2015also Dennison and Ogilvie 2014).
Altogether, our results forcefully suggest that historical Europe had several systems of household formation, and not just two. Whereas our clusters 1 and 4 (k ¼ 4) encapsulate well the features Hajnal had in mind when contrasting his two systems of household formation, the two other partitions that we have accounted for are impossible to classify within the Hajnalian framework, because the characteristics of these groupings cannot be limited to a gradient in a gradation line such as in an intermediate "transitional" classification (see Laslett 1983;cf. Micheli 2018). In fact, our cluster-based taxonomy should be viewed as competing configurations of the four interlinked variables, rather than as the gradual transition from one extremity to the other.
The spatial spread of the four-structure partition presents a much more complex picture than any of the "lines" or "zones" that Hajnal and other researchers have previously been able to map out. Neither the "unique" household formation type that Hajnal was talking about, nor its opposite (clusters 4 and 1 mentioned above, respectively), were found to be fully consistent with the geographic areas for which their respective "ideal types" were said to be a norm. The most extreme version of the Hajnalian "western" household formation pattern was detected in part of Scandinavia and northern-central Germany; but not in England, the Netherlands, or northern France, where it might be expected to prevail (Hajnal 1982, 449; for similar observations, see Dennison and Ogilvie 2014). What is more, "western-like" household formation rules also prevailed among indigenous peoples of the Russian north (Siberia; the shores of the Barents Sea), that is well beyond the supposed line dividing "European" and "non-European" family regimes (Hajnal 1965(Hajnal , 1982. The "joint family" cluster is spatially discontinuous as well, with its eastern European and Balkan areas being separated from one another by territories where household formation was based on entirely different principles. The spatial distribution of Cluster 2 also shows that the latter's attributes are not confined to any particular "cultural" region of Europe (cf. Macfarlane 1981), while the geography of Cluster 3 makes it clear that the characteristics that have long been embraced as the historical metric of north-western European "exceptionalism" (Hajnal 1982;Laslett 1977;Smith 1993; de Moor and van Zanden 2010) may have been shared with many other European populations (Goody 1996;Ruggles 2009;cf. Hajnal 1982, 450).
Taken together, the four-cluster structure established in this paper offers two main lessons for further comparative family history research. By showing that the structural progression within the areas covered by our datai.e., away from more nucleated, and more neolocal households with high incidences of service and late marriage; and toward populations with earlier marriage, weaker marriage-headship nexus, more household complexity, and less servicedoes not consistently travel along the west-east axis, our results suggest that the habitual ways of conceptualizing the European familial past in strictly bipolar categories need to be abandoned. By the same token, it also warns against too hasty comparisons between Europe and Asia (cf. Lundh and Kurosu 2014). As no standardized "dual Europe" can be discerned from the historical data we have analyzed, any gross comparisons between the situations across the Eurasian landmass must be approached with caution (cf. Goody 1996).
Finally, this study has important limitations that lay the ground for subsequent research. Since our results are first of the kind to appear, unresolved remains the question of what might lie behind the revealed differences in household patterns in different parts of Europe. Although an exploration of these issues is beyond the scope of this paper, it is clear that the four-cluster structure extends beyond potential cultural, national, or even broad regional boundaries, thus calling for complex explanations that take into account a combination of historical, political-economic, and ecological factors (Kertzer 1991a;Rudolph 1992). Future comparative research will have to take a stance toward this nuanced geography and use appropriate tools to unravel the nature and permeability of "frontiers" or "borderlines" in factors underlying household formation systems in the past.
Second, in this paper we have no data on Italy and the Iberian peninsula (except for Catalonia); that is, for regions that would likely be found to display a wide range of household formation patterns (Barbagli 1991;Rowland 2002). However, we have good reasons to expect that filling this gapif doing so is ever possiblewould only strengthen our results, as at least the Italian sub-regional variation in household formation patterns would not reveal a specific distinct monolith-like cluster within our data structure, but would instead be partitioned between the existing cluster structure that we have derived. Nevertheless, further research including Italy might prove useful, given that part of its suggested micro-demographic variability has recently been claimed to have been forged on the back of geographical biases in data selection (Curtis 2015). As for now, however, our results confirm that the claim made about Italy-i.e., that it is "a burial ground for many of the most ambitious and well-known theories of household and marriage systems" (Kertzer 1991b, 247) may have wider pan-European implications, and cannot be dismissed as a mere peculiarity of the Mediterranean.
Another potential criticism of our work might be that our clustering algorithms searched for similarities and dissimilarities between populations from different points in time. Like the rich historiographical tradition our approach was based on, 24 we faced some obvious challenges: namely, that for the reconstruction of historical household formation patterns on a large pan-European scale, an ideal data structure that includes cross-sectional data from around the same point in time, or continuous time-series for all populations in our database does not yet exist, and is unlikely to become available in the future. As for the current data structure, it has to be pointed out that although the distribution of the four (base) clusters between time periods is not equal, most of these clusters contain regional populations that cover the three time periods that our data span. This may suggest that we should be careful in assuming that the partitions we have identified are driven by data conglomerates that are strictly time-bound. Although we cannot rule out the possibility that a different chronological composition of our dataset would have produced different partitions, the ability of such data permutations to seriously undermine our findings remains doubtful, at least within the limits of currently available data. 25 Still, our reliance on pooled-time cross-sections may rise problems in so far as it concerns the large samples drawn from the 1851 census for Britain (and the 1881 census for Scotland). The rural communities in these censuses were economically far more differentiated than in most of Europe and the farming methods most of them relied upon were more capitalistic then elsewhere (e.g., Overton 1996;Gritt 2002), principally employing wage labor with farm servants long a feature of the past. Had there been data for Britain 100-150 years earlier, it might be argued that the servant component could have been more prominent (and marriage ages could have been very different). 26 However, to the best of our knowledge, such data are not available in a computerized form (Wall, Woollard, and Moring 2004), but even if they were, they would miss information crucial in the computation of the household formation markers (ages in particular). Moreover, the extent to which the inclusion of these data would alter the cluster-based structure presented above, particularly by pushing England apart from east-central European populations (e.g., Poland) remains uncertain. Szołtysek (2015), for example, has shown that the percentage of servants and the percentage of households employing them in 18th-century Poland-Lithuania were, respectively, almost identical or significantly higher than the comparable shares in England as captured with Laslett's "master" sample of the 63 early modern parishes (Szołtysek 2015, 337).
We should also not ignore the hermeneutic limits of our measures and of the taxonomies that they allowed us to derive. One of the well-recognized limitations of Hajnal's household formation hypothesis was that he did not include stem families into his typology (see Engelen and Wolf 2005, 16-18;Szołtysek 2015, 574-76). It might be argued that our following Hajnal in this respect pose a threat to the validity of our groupings by suggesting that in some areas important variations within non-nuclear family households may be missed. Nevertheless, we embrace the view that for household formation theory to be improved in the future, information complying with Hajnal's theoretical background ought to be used in the first place to produce comparative research despite its flaws. It should also be remembered that there is a widespread disagreement about exactly how to define the stem-family system and that pursuing this task with quantitative methods alone is predestined to fail because meaningful indications of the presence or absence of the stem-family pattern should take into account behavior and norms, as well as the notion of sequence (Szołtysek 2016). Naturally, these difficulties do not exempt us from seeking new ways of capturing the stem family in a broader comparative perspective, but this task must be left for further research (cf. Szołtysek et al. 2019). The results presented above offer a useful and necessary starting point for such explorations and we hope they will provide a stimulus for inquiring about where and how the stem-family configurations could be accommodated into a reliable taxonomy of European household formation systems.
Finally, it is obvious that Hajnal's four household formation markers do not cover all of the axial principles of the family systems in historic Europe, and even less so beyond it. There are other elements that may need to be taken into account when investigating the differences and the commonalities among the local family systems within Europe (e.g., Gruber and Szołtysek 2016). Future research should strive to more directly establish the affinity of the systemic and spatial relationships presented in this paper with those more complex classificatory ventures. The intellectual challenge of meaningfully assessing the variability of historical family patterns in Europe is still there.

Notes
1. An abridged version of the paper was published under the same title in Wall and Robin (1983). For the sake of conceptual clarity, this paper sticks as closely as possible to Hajnal' Carmichael et al. 2016). 2. In philosophy, the term "natural kinds" refers to the idea that there are some "naturally" separated classes in observer-independent reality, which, for traditional realists, may correspond to "true clusters" (Hennig 2015, 54). In mathematical formalism, "naturalness" implies that members of natural groups are more closely related to fellow members than they are to non-members (Jain 2010 (1982,452) were presented as applying to both sexes, in this paper, marriage and service variables are computed for women only. Male and female variables are highly correlated across our data (Pearson's r¼ .806, significant at the 0.01 level for marriage; and, r ¼ .832; significant at the 0.01 level, for service, accordingly). However, Hajnal himself, as well as many of his followers, paid more attention to the variability of marital ages among women rather than among men (Hajnal 1965, 134). Furthermore, the incidence of female servants is a more sensitive gauge of social admissibility of service, given the historically generalised constraints on female autonomy (Kok 2017;Gruber and Szołtysek 2016). Following Hajnal (1982, 463-466), the headship attainment has been computed for males only. 8. "Servants" were identified using the variable "relationship to head of household," which indicated that these individuals were live-in servants, apprentices, or journeymen (or individuals whose designations can be translated with these words). 9. This is an option of the DescTools package (see Signorell et al. 2018). 10. Given that clustering algorithms can find clusters even if no meaningful groups are embedded in the data (Ketchen and Shook 1996), our dataset was first tested for the presence of a uniform (random) data distribution. The Hopkins statistic (Lawson and Jurs 1990), which yielded a value close to zero (0.1612183), suggested that our dataset consists of significantly clusterable data (Kassambara 2017, 123-124). The visual assessment of the cluster tendency (VAT; Kassambara 2017, 124; figure available upon request) also confirms that there is a cluster structure in the data set. Since the variables we selected are measured on a different scale, each of the four metric variables were z-standardised prior to the analysis for this and for all subsequent analyses. 11. The proceeding cluster analysis was carried out using the R package 'cluster' (Maechler 2018). 12. The basic idea of this algorithm is to first compute randomly selected medoids from the K y data objects to form the K y cluster. After finding the set of medoids, each object of the dataset is assigned to the nearest medoid using the Euclidean distance metric. That is, object i is put into cluster v i when medoid mv i is nearer than any other medoid mw. The whole process is repeated to ensure that the objective function (E; the sum-of-squared-error between the objects in the clusters and their respective cluster centroids) cannot be reduced by interchanging (swapping) a selected object with an unselected object. This procedure is carried out iteratively until none of the data points move to a different cluster (local minimum). At this stage, the clustering algorithm is said to be complete (Kaufman and Rousseeuw 1990;Kassambara 2017, 48-56;Han et al. 2012, 454-457. 13. The latter method is also inferior to PAM in terms of sensitivity to outliers. 14. NbClust R package was used (Charrad et al. 2014). 15. Note, that for both types of partitions, the estimated silhouette width (0.57 and 0.36 respectively) represent either "good" or "fair" clustering solution (Mooi and Sarstedt 2011, 280;also Finch et al. 2015). 16. Altogether, 85.3 per cent of the original information has been retained in two-dimensional space. Thus, the obtained dimensionality reduction implies very little information loss. 17. The first and the second cluster in this partition each covers 12 per cent of all of the regions. The fourth cluster covers around 22 per cent of the regions, while the third cluster covers 54 per cent of the regions. 18. Both the spatial and the temporal dispersion of this cluster membership suggest that its geography cannot be attributed solely to either the census chronology or the temporal trajectory of the industrial revolution. 19. The dimensionality reduction reveals that almost 86 per cent of the point variability is explained by the first two principal components. 20. It has been documented that 65 percent of domestic groups in the nuclear mode may well indicate the presence of the stem-family rules under pretransitional demographic conditions (see Berkner 1976). In the stem family mode, only one of the children marries and lives with his/her parents in a multigenerational household, while his/her siblings have to either leave or stay unmarried while remaining in the parental household. If such a custom additionally implies the takeover of the farm by the young married couple, the CMHD would still be rather low, whereas the overall share of nuclear families would decrease (see Hajnal 1982, 453). 21. The risk of this being caused by the predominance of post-1850 populations in the cluster must be ruled out. Our service variable primarily captures the phenomenon of women-driven domestic service, the prevalence of which wasto the best of our knowledgequite robust to the impact of 19thcentury modernisation processes, as it either remained high or declined very little in this period (Hinde 1985;also Woodward 2000;Gritt 2000). See, however, also the Discussion and Conclusion section. 22. We perform the analysis with the use of R package 'clValid', see: Brock et al. 2008. Only APN is normalised in the interval (0.1); the remaining indices need to be minimised. 23. The method is based on a comparison between every single cluster in a clustering and its "mirror image", which is obtained with the use of four resampling methods that are compared by means of the Monte Carlo simulation with 1,000 resampling runs (Hennig 2007, 258). 24. Laslett's well-known "English sample of 64 settlements" used widely for the comparative study of domestic groups has dated from 1574 to 1821. In a similar vein, Hajnal (1982, 455) used Klapisch/ Herlihy's data on the medieval Florence catasto to argue that 15th-century Tuscany displayed a joint family formation system that was remarkably similar to those documented for 20th-century China and India. In addition, to illuminate the differences in the marriage-headship nexus between different family systems, Hajnal compared data from Denmark in 1801, and from the Pisa area in 1427-30. Dennison-Ogilvie's recently created meta-database (Dennison and Ogilvie 2014) also includes data from different time periods, and these data are used to compile country or regional averages. 25. A number of tests have been performed in which the NAPP input data were changed so that they could, in some cases, replace the oldest censuses used in the base clustering with later censuses (i.e., Iceland 1701 was changed into 1801; Denmark 1787 was changed into 1801; Norway 1801 was changed into 1865; England and Wales 1851 was changed into 1881, and Sweden 1880 was changed into 1890). Despite such aggressive changes being introduced (i.e., "traditionality" or "rurality" are much harder to assume in this case than they are for the base input data structure), the overall level of agreement of the two methods of clustering is surprisingly high: for k ¼ 2, an almost perfect overlap with the base partition was found (AMI ¼ 0.99); for k ¼ 4, a moderate overlap for the two methods (AMI ¼ 0.50) was observed. 26. The long-term persistence of the major features of the English domestic group (in particular it's nuclear structure) has been proclaimed by a number of prominent scholars since the publication of Laslett's seminal overview (Laslett 1969;e.g., Wall 1977). Major residential patterns in England and Wales changed relatively little also over the period 1851-1911(Sch€ urer et al. 2018).